      subroutine blufac

      use param

#include <global.h>
      include 'ntbytes.h'     
      include 'geom.h'
      include 'setls.h'
      include 'cmfd.h'

      real(NBF) rdet 

      do im=1,ng
        bdd(:,1,im)=diagf(:,1,im)
        bdu(:,1,im)=cczf(:,2,1,im)
        do j=2,tnode
          ja=j-1
          rdet=1/(bdd(1,ja,im)*bdd(4,ja,im)-bdd(2,ja,im)*bdd(3,ja,im))
!     set lower block
          bdl(1,j,im)=rdet*(cczf(1,1,j,im)*bdd(4,ja,im)- &
          cczf(2,1,j,im)*bdd(3,ja,im))
          bdl(2,j,im)=rdet*(cczf(2,1,j,im)*bdd(1,ja,im)-   &
          cczf(1,1,j,im)*bdd(2,ja,im))
         bdl(3,j,im)=rdet*(cczf(3,1,j,im)*bdd(4,ja,im)- &
                     cczf(4,1,j,im)*bdd(3,ja,im))
          bdl(4,j,im)=rdet*(cczf(4,1,j,im)*bdd(1,ja,im)- &
                     cczf(3,1,j,im)*bdd(2,ja,im))
!     set diagonal block
          bdd(1,j,im)=diagf(1,j,im)-bdl(1,j,im)*bdu(1,ja,im)- &
                     bdl(2,j,im)*bdu(3,ja,im)
          bdd(2,j,im)=diagf(2,j,im)-bdl(1,j,im)*bdu(2,ja,im)- &
                     bdl(2,j,im)*bdu(4,ja,im)
          bdd(3,j,im)=diagf(3,j,im)-bdl(3,j,im)*bdu(1,ja,im)- &
                     bdl(4,j,im)*bdu(3,ja,im)
          bdd(4,j,im)=diagf(4,j,im)-bdl(3,j,im)*bdu(2,ja,im)- &
                     bdl(4,j,im)*bdu(4,ja,im)
!     set upper block
          bdu(1:4,j,im)=cczf(1:4,2,j,im)
        enddo
      enddo
      
      end subroutine
